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A gap in understanding the link between continuum theories of ion transport in ionic liquids and 
the underlying microscopic dynamics has hindered the development of frameworks for transport 
phenomena in these concentrated electrolytes. Here, we construct a continuum theory for ion 
transport in ionic liquids by coarse graining a simple exclusion process of interacting particles on a 
lattice. The resulting dynamical equations can be written as a gradient flow with a mobility matrix 
that vanishes at high densities. This form of the mobility matrix gives rise to a charging behaviour 
that is different to the one known for electrolytic solutions, but which agrees qualitatively with the 
phenomenology observed in experiments and simulations. 

PACS numbers: 82.45.Gj, 82.45.Fk, 47.57.jd, 78.30.cd 


Room temperature ionic liquids play an increasingly 
important role as electrolytes in electrochemical and elec¬ 
tromechanical applications ranging from actuators [lj-@] 
to supercapacitors [JI3- Ionic liquids differ from tradi¬ 
tional electrolytes in that they consist only of positive 
and negative ions without any solvent. Recent theoret¬ 
ical calculations jsj suggest that ionic liquids are con¬ 
centrated electrolytes, and should not be modelled as a 
weak electrolyte (with solvent consisting of ion-pairs), as 

suggested elsewhere 113- 

In many important technological applications, ionic 
liquids are close to an electrified interface [11]. Although 
the equilibrium structure of the electrical double layer 
is relatively well studied, understanding the dynamic re¬ 
sponse of ionic liquids to an applied potential or sur¬ 
face charge is more challenging because of the difficulty 
in identifying an appropriate non-equilibrium dynamic 
framework. Previous theoretical studies 12|- 17] relied on 
the dynamical density functional theory 118l. 1 191]. in which 
the ion flux, j±, is related to the ion density fields c± and 
free energy density functional ,F[c±] via 


„, „ , SF 

J± = —M±c±V [ — 


(1) 


where M± is the cation/anion mobility. A key as¬ 
sumption in the derivation of © is that ion density is 
low compared to an underlying solvent bath [_2Cj]. This 
assumption may become problematic, as can be illus¬ 
trated by substituting the lattice gas free energy, T ex = 
(k B T/a 3 ) / f2 [a 3 clog(a 3 c) + (1 - a 3 c) log(l - a 3 c)] d 3 z 
where fl is system’s volume, into ©; applying the con¬ 
tinuity equation, we obtain 


dc 

dt 


= DV • 


Vc 


1 — a 3 c 


( 2 ) 


where D = M/k B T is the diffusion constant, k B the 
Boltzmann constant and T temperature. However, the 
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FIG. 1: Schematic of the system under consideration: cations 
and anions on a lattice of lattice constant a. 


continuum limit of a system of particles on a lattice with 
lattice constant a undergoing a simple exclusion process 
is well known, and leads to the linear diffusion equa¬ 
tion dc/dt = DS7 2 c for the particle density c 21], rather 
than ©. 

Lattice gas models of ions, first proposed by Bikerman 
f22| in the 1940s, are commonly used as simple models 
of ionic liquids 12, [3, liij - Ebl ] in equilibrium. The goal 
of this paper is to derive a consistent model for the dy¬ 
namics of ions in solvent-free ionic liquids, and analyse 
the dynamics of electrical double layer formation. We 
map the system onto a lattice and take the continuum 
limit of the microscopic reference kinetics of a discrete 
symmetric exclusion process (see Figure©. This refer¬ 
ence kinetics is a natural one to consider as ion motion 
in a concentrated assembly is physically akin to particles 
“hopping” on a lattice 45[. Similar reference kinetics 
were successfully used to model spinodal decomposition 
in alloys and were shown to be a microscopic ba¬ 

sis for the Cahn-Hilliard equation [IM1. Our approach 
has the advantage that steric exclusion is accounted for 
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at the level of dynamics. 

We first consider a one-dimensional lattice of lattice 
constant a (corresponding to ion diameter a in the contin¬ 
uum limit) for simplicity, and later generalize our results 
to higher dimensions. We consider a discrete-time dy¬ 
namics in which particles can only move between nearest- 
neighbour lattice sites between time t and t + At. De¬ 
noting by S“(i) e {0,1} the occupancy of the i th lattice 
site at time t by ion of type a = {+,—}, the evolution 
master equation for reads 

ST(t + At) = r^ i+1 s?s i+1 + 7? +1 _*S? +1 ( 1 - S t ) 

+ rfLu^SfSi -1 + - 5 <) 

+ (1 - rUi+i - rUi- 1) Sf, (3) 

where the Sf on the right hand side are taken at time t 
and Si = Sf + S~\ r^ i± i is unity if particle a at site i 
attempts to jump to site i ± 1 and is zero otherwise: this 
transition propensity takes into account the long-ranged 
interparticle interactions as explained below. The first 
term of Eq. © ensures that there will be particle a at 
site i and time t + At, if it is there at time t, attempts 
to move to site i + 1 and finds that site occupied. The 
second term describes a possible transition from site * +1 
to site i when there is no particle at site i. The third and 
forth terms describe the same processes between sites i 
and i— 1. Finally, the last term corresponds to a particle 
a at site i that does not attempt to leave it during the 
interval At. 

An ensemble average of with \i — j\ = 1 gives a 
transition rate that satisfies the detailed balance and can 
thus be related to the Boltzmann factor by 

(.£*> = i (4 ) 

where 

r = (5) 

P 

is a potential acting on particle a at site i due to 
all remaining particles, and U a p(\i — j\) is the micro¬ 
scopic (electrostatic) interaction potential between the 
particles. In the absence of long-ranged interactions 
(U a p = 0), there is no preferred direction of motion, thus 

= V 2 - 

The continuum evolution equation can be obtained 
by rescaling the lattice indices by the lattice spacing a, 
and introducing the minimal lattice volume v (v = a 
in 1-D) as well as the ensemble average concentrations 
c a (x = ai) = v^ 1 (Sf) and mean potentials p l a (x) = 

3 fn U a p(\x — x'\)cp(x')dx'. Taking the average of both 
sides of Eq. ©, and applying the mean field approxima¬ 
tion, (S!j*Sj • • •) « (SP‘)(Sj) • • •, we can expand the re¬ 


sulting expression in a power series in a and At to obtain 


1 dc a d 2 c d 2 c a 1 d 2 n l a 

= vc nTr -;-{-(l-vc)— T +c a {l-vc)- 


D dt 


‘ dx 2 


dx 2 


dc dc a 

VC -^Z ~ i 1 ~ vc >^~ 


dx 


dx 


k B T dx 2 
1 dn\ 


k B T dx 


( 6 ) 


where we have defined c = c+ + c_ as the total ion den¬ 
sity, and identified D = a 2 / (2At) as the self-diffusion 
coefficient. 

To generalize Eq. © to higher dimensions, we assume 
that the fluxes along different axes are decoupled. Intro¬ 
ducing the mobility matrix 


M = 


D ( c- t-(l — vc ) 


0 


k B T \ 0 c_(l — vc) 

the higher-dimensional version of © is 
dc 


dt 


= V-CMV/i), 


( 7 ) 


( 8 ) 


where c T = (c+,c_), = (/x+,/x_), and is defined 

by = 8J : /5c± where 

F[c±] = 7T / c a (xi)[/ a/3 (|xi - x 2 |)c / 3(x 2 ) clxidx 2 
z . q J n 


k B T r 
v Jn 


vc a log(wc a ) + (1 — vc) ln(l — vc) 


dx. 


(9) 


Equation © is the continuum kinetic equation for an 
interacting two component system. Note the important 
physical constraint that the evolution equation for the to¬ 
tal concentration, c, in the absence of long-ranged inter¬ 
actions (U a p = 0), reduces to the linear diffusion equa¬ 
tion. This constraint, as explained above, respects the 
fact that the underlying dynamics of our reference sys¬ 
tem is a simple exclusion process on a lattice. Contin¬ 
uum kinetic equations with the same mobility function as 
© have been proposed in the literature in the context of 
the modified Cahn-Hilliard eqi 
models of Li-ion batteries 
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uation 


34J, and phase-field 


To apply Equation © to an ionic liquid system, we 
introduce a characteristic length scale l c for short-ranged 
interactions, and split the Coulomb potential U a p = 
U S J R + U l J 0 , where U^{x) = q a qpl B e~ x / 1 -/x, U lx p (x) = 
_ e -x/i c 


)/x |37|,[38|, and l B = e 2 /(47reo ek B T) 


J ap _r 

q a qplB{X 

is the Bjerrum length. Below the length scale l c , it is ac¬ 
tually the hard core exclusion that matters rather than 
Coulomb interaction, thus U^ p can be neglected. This 
truncation of the Coulomb potential is necessary as our 
mean-field approach underestimates steric correlations, 
and as such the divergence of the Coulomb interaction 
at the origin renders electrostatic interactions effectively 
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too strong. We thus write 


U aP {x) U l ^{x) 1 — e~ x / lc 

~kWT~ ~ ~k^T~ = QaqplB - x - 


( 10 ) 


This decomposition of the Coulomb potential is not 
unique — the exponential function is chosen phenomeno¬ 
logically and for mathematical convenience. 

Introducing the local electric field u , and exploiting the 
Green function, one can rewrite the non-local integro- 
differential equation (j8]) as a set of coupled partial differ¬ 
ential equations 


(1-^V 2 )V 2 m = 


dc± 

~dt 


DX7 • c±(l 


—47 tIb(c+ — c_), 


— vc)V 


±.u + In 


vc± 

1 — VC 


( 11 ) 

( 12 ) 


Equation CD is identical to the modified Poisson equa¬ 
tion derived phenomenologically in [ 2 f| using a gradient 
expansion of a nonlocal electrostatic kernel. We note 
that [Dot j took the variational approach of [Hi to develop 
a framework for charge-transfer reaction kinetics, with 
the resulting equation similar to (fl2l) . Here we provided 
a microscopic statistical derivation of the kinetics of ion 
transport. 

We turn our attention to a simple problem to gain 
some insight into the characteristic behaviour of CD- 
©: an ionic liquid with bulk cation/anion concentra¬ 
tion Co bounded by two parallel, blocking electrodes at 
x = —L, L. Initially the concentrations of the two ion 
species are uniform, and a step voltage of amplitude 2V 
is applied at t = 0 + . Introducing the Debye length 
Id = l/y/Ei tcqTb and dimensionless packing parameter 
in the bulk |12L |23( 7 = vco, we introduce the dimen¬ 
sionless variables r = ( D/Llo)t , X = x/L , C± = c±/cq. 
The no-flux conditions at the electrodes read 


,, ^ du , „.dC± „ dC 

±c±(1 " lC) dx + (1 - 7C) ^r + lC± ox 


= 0 . 


J x=±i 
(13) 

At the electrodes surface, we posit that the classical 
Gauss law ±eu x = 47 rer holds at X = ± 1 , with a the 
(dimensional) surface c harg e density, and e the dielectric 
constant of the medium mm. This condition, together 
with the constant potential condition gives 


u{X = ±l,r) = ±V, u X xx(X = ±1,t) = 0, r > 0. 

(14) 

The initial conditions are 


c±(i,o) = i,ie[-i,i], (15) 

To avoid complications of double layer overlap, we con¬ 
sider widely separated electrodes taking L/Id = 100 . We 
take l c = a = i? 1 / 3 , Bjerrum length Ib = 50A, ion diam¬ 
eter a = 5A and 7 = 0.25 (see e.g. [HI though the 
qualitative behaviour reported below is not sensitive to 
7 ); we therefore have I c /Id = ^8777 Ib/ a « 10 . 
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FIG. 2: (a) The total density, C++C-, and (b) charge density 
C+ — C-, as functions of distance X = x/L from the elec¬ 
trode and time after an applied step voltage V = 40Vr with 
Vt = fcsT/e. Here 7 = 0.25, I c /Id = 10 , and L/Id = 100 , 
red/blue curves denote the first/second charging regimes, and 
the arrow shows the increasing density deficit in the first 
charging regime. Numerical solution of Eqs. CED-CE! is per¬ 
formed using pdepe in Matlab. 


Figure [2] shows that charging proceeds through two 
distinct regimes: First, the (negative) electrode attracts 
cations from the vicinity and expels anions, resulting in a 
dense, “compact layer” of cations near the electrode that 
overcompensates the surface charge (region I in Figure 
n. Ion diffusion is hindered as the mobility matrix 0 
vanishes in regions of high density. As a result, the total 
density reaches a minimum away from the compact layer 
(c./. red arrow in Figure[2^)- In the second stage, anions 
arrive from the bulk to screen the now net-positive com¬ 
pact layer. This flux fills the total density deficit near 
the compact layer incurred in the first charging regime, 
creating a region of negative charge density and in fact 
excess total density (region II in Figure [DO- 

A key measure of practical interest is the integrated 
total diffuse charge, 

Q(r) = J \C+(X,t) - C-(X,t)] dX. (16) 

Note that the overall system is electroneutral, therefore 
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FIG. 3: (Main figure) The total charge as a function of time 
for different values of I c /Id with L/Id = 100 , 7 = 0.25 fixed. 
(Inset) The equilibrium relaxation time r re i ax as a function of 
l c and Id (obtained by fitting numerical results to an expo¬ 
nential decay). 



FIG. 4: The voltage drop across the system evolves non- 
monotonically under constant current conditions, 05). Here 
L/Id = 100 , Ic/Id = 10 and 7 = 0.25, and the total 
charge Q = Jr. Inset: simulation data from [42] for which 
« 400 kAcnG 2 corresponds to dimensionless J = 2. 


the total charge of the ions is equal and opposite to the 
surface charge. Q is therefore the charge accumulated 
at the anode, which is equal and opposite in sign to 
the charge accumulated at the cathode. Figure [3] shows 
that, as charging proceeds, the total charge initially in¬ 
creases, corresponding to the formation of the compact 
layer. However, arrival of anions in the second charg¬ 
ing regime decreases the charge to the final equilibrium 
value. This charging mechanism is schematically illus¬ 
trated in the inset of Figure [3] The correlation length 
l c /1 d controls the extent of charge oscillation and thus of 
overcompensation of electrode surface charge by the com¬ 
pact layer. Therefore, decreasing the correlation length 
reduces the extent of charge overcompensation and also 
the peak diffuse charge. 

Further insights into the charging process can be ob¬ 
tained by noting that the initial rise in charge occurs 
over r = 0(1). In dimensional terms this cor resp onds 
to t,Rc = LId/D , the usual RC time constant 4l|, cor¬ 
roborating the fact that the peak has its origin in the 
formation of the diffuse layer. Numerical experimenta¬ 
tion (see inset of Figure [3| suggests that the late-stage 
exponential relaxation of the charge to equilibrium has a 
distinctly different timescale 


irelax ~ oil. 


3/2 


(17) 


This scaling suggests that the decay in the stored charge 
comes from the formation of charge oscillations: L 2 /D 
gives the decay time due to diffusion of ions through the 
electrochemical cell, and this is rescaled by {Id/Ic) 3 ^ 2 
where k c ~ is the characteristic wavelength of charge 
oscillations (c.f. Equation (HD). 

The non-montonic evolution of Q(t) is in stark contrast 
to the results predicted by dynamical density functional 


theory [3IISE3, where the diffuse charge is mono- 
tonically increasing. We note that this effect is different 
from kinetic charge inversion due to double layer overlap 
[l5l ] . The degenerate mobility ([7|) in our approach en¬ 
sures that the flux due to electrostatic interactions van¬ 
ishes at close packing, and thus there are distinct regimes 
of initial charge density polarisation and, at later times, 
rearrangement of the double layer into cation-rich and 
anion-rich layers. 

Qualitatively similar behaviour is obtained under 
charge-controlled conditions, i.e. imposing a constant 
current, 


du 

dX 


= ± Jr. 


(18) 


-Y=±l 


Figure [4] shows that the non-equilibrium double layer re¬ 
arrangement manifests itself in the non-monotonic evo¬ 
lution of the potential drop across the system when the 
current density J is large. This qualitatively agrees with 
recent molecular dynamics simulations 42], but is in con¬ 
trast to conventional dynamical density functional the¬ 
ory, which again predicts a monotonic increase in poten¬ 
tial drop as a function of time. 

In summary, we have derived a continuous model for 
the dynamics of solvent-free ionic liquids based on coarse- 
graining a simple exclusion process of interacting parti¬ 
cles defined on a lattice. The resulting equations have 
the structure of a gradient flow with a degenerate mobil¬ 
ity function. As examples, these equations were analysed 
for a system where: (i) a step voltage is applied between 
widely separated electrodes, and (ii) a constant charging 
current is applied. Even in these simple cases, our the¬ 
ory differs qualitatively from previously developed the¬ 
ories for electrolyte solutions. Importantly, we showed 
that the total diffuse charge is a non-monotonic function 
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of time. Experiments and simulations of the dynamics 
of ion transport in ionic liquids are currently scarce; we 
hope that our theory provides a framework to interpret 
experiments and motivate further investigation. 
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